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Abstract 

This paper presents a theoretical discussion as well as novel solution algorithms 
for problems of scattering on smooth two-dimensional domains under Zaremba bound¬ 
ary conditions—for which Dirichlet and Neumann conditions are specified on various 
portions of the domain boundary. The theoretical basis of the proposed numerical 
methods, which is provided for the first time in the present contribution, concerns de¬ 
tailed information about the singularity structure of solutions of the Helmholtz operator 
under boundary conditions of Zaremba type. The new numerical method is based on 
use of Green functions and integral equations, and it relies on the Fourier Continua¬ 
tion method for regularization of all smooth-domain Zaremba singularities as well as 
newly derived quadrature rules which give rise to high-order convergence even around 
Zaremba singular points. As demonstrated in this paper, the resulting algorithms enjoy 
high-order convergence, and they can be used to efficiently solve challenging Helmholtz 
boundary value problems and Laplace eigenvalue problems with high-order accuracy. 
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1 Introduction 


This paper concerns problems of scattering on smooth two-dimensional domains under 
Zaremba boundary conditions, and, in particular, it presents a theoretical discussion as 
well as novel solution algorithms for this problem. The theoretical basis of the proposed 
numerical methods concerns detailed information, put forth in Section 4, about the singu¬ 
larity structure of the solutions of the Helmholtz operator under boundary conditions of 
Zaremba type. The new numerical method, which is based on use of Green functions and 
integral equations, incorporates use of the Fourier Continuation method for regularization 
of all smooth-domain Zaremba singularities as well as newly derived quadrature rules which 
give rise to high-order convergence even around Zaremba singular points. As demonstrated 
in this paper, the resulting algorithms enjoy high-order convergence, and they can used 
to efficiently solve challenging Helmholtz boundary value problems and Laplace eigenvalue 
problems with high-order accuracy. 

We thus consider the classical Zaremba boundary value problem 


Au{x) + k‘^u{x) 
u{x) 
du{x) 
drix 


0 

X G D, 

f{x) 

X G Fj) 

9{x) 

X G Fjv 


( 1 . 1 ) 


for the Helmholtz equation, where H is a 2-dimensional domain with boundary T consisting 
of two disjoint portions T o and Tjv upon which Dirichlet and Neumann values are prescribed, 
respectively; see Figure 1. 
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Figure 1: Generic smooth domain G with boundary partitioned into Dirichlet and Neumann 
regions shown in red and blue, respectively. 


Signihcant challenges arise in connection with Zaremba boundary value problems in view 
of the singular character of its solutions: as hrst established by Fichera [14], Zaremba solu¬ 
tions are generally non-smooth even for infinitely differentiable boundary data / and g, and 
smoothness of solutions can only be ensured provided / and g obey certain relations which 
generically are not satisfied. Wigley [29, 30] provides a detailed description of the singularity 
structure around Zaremba points—which includes singularities of both square-root and log¬ 
arithmic type. In Section 4 of the present paper it is shown, however, that a tighter result 
holds in the case the domain boundary is itself smooth. 
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Problems with mixed Dirichlet-Neumann boundary conditions were first considered in 
Zaremba’s 1910 contribution [32], which established existence and uniqueness of solution in 
the particular case of the Laplace equation (/c = 0). Boundary conditions of Zaremba type 
arise in a number of important areas, including elasticity theory (were it appears as a model 
in the contexts of contact mechanics [28] and crack theory [12]); homogenization theory (as it 
applies to problems of steady state diffusion through perforated membranes [12]), etc. One of 
the main motivations for our consideration of this problem concerns computational electro¬ 
magnetism: the Zaremba problem serves as a valuable stepping stone to the closely related 
but more complex problem of electromagnetic propagation and scattering at and around 
structures such as printed circuit-boards. (As in the Zaremba problem, where the boundary 
conditions change type at the Dirichlet-Neumann boundary, the boundary conditions for the 
Maxwell equations at and around circuit-boards change type—not from Dirichlet to Neu¬ 
mann but from dielectric transmission conditions to perfect-conductor conditions—at the 
edges of the perfectly condncting circnit elements.) 

After the initial contribntion by Zaremba, early works concerning the Zaremba bonnd- 
ary valne problem inclnde results by Signorini [25] (1916: solution of the Zaremba problem 
in the npper half plane using complex variable methods); Giraud [17] (1934: existence of 
solntion of Zaremba problems for general elliptic operators); Fichera [13,14] (1949, 1952: 
regularity studies at Zaremba points, Zaremba-type problem for the elasticity equations in 
two spatial dimensions); Magenes [22] (1955: proof of existence and uniqneness, single layer 
potential representation); Lorenzi [21] (1975: Sobolev regularity around a corner which is 
also a Dirichlet-Nenmann junction); and Wigley [29,30] (1964, 1970: explicit asymptotic 
expansions aronnd Dirichlet-Nenmann junctions), amongst others. More recent contribn- 
tions in this area include reference [28], which provides a valuable review in addition to a 
stndy of Zaremba singnlarities and theoretical resnlts concerning Galerkin-based computa¬ 
tional approaches; reference [7], which considers the Zaremba problem for the biharmonic 
eqnation; references [10,11], which study Zaremba boundary valne problems for Helmholtz 
and Laplace-Beltrami eqnations; reference [8], which discusses the solvability of the Zaremba 
problem from the point of view of pseudo-differential calculus and Sobolev regnlarity the¬ 
ory; reference [18] which introduces a certain inverse preconditioning technique to reduce 
the number of linear algebra iterations for the iterative nnmerical solution of this problem 
and which gives rise to high-order convergence; and hnally, reference [4], which successfully 
applies the method of difference potentials to the variable-coefficient Zaremba problem, with 
convergence order approximately equal to one. 

The rest of this paper is organized as follows: after some necessary preliminaries presented 
in Section 2, Section 3 lays down the integral eqnation system we use, and it studies the 
equivalence between the integral eqnation system and the original Zaremba problem (1.1). 
Detailed asymptotics for the Zaremba singularities at Dirichlet-Neumann junction are pre¬ 
sented in Section 4. Bnilding upon these results, further, and exploiting an interesting 
connection of this problem with the Fonrier-Gontinuation method [3,6], Fonrier expansions 
for the integral densities are obtained in Section 5 which regularize all Zaremba singulari¬ 
ties. Section 6 then introduces a numerical algorithm that incorporates the aforementioned 
Fourier Gontinuation series on the basis of certain closed form integral expressions, and 
which thus produces numerical solutions with high-order accuracy. A variety of results in 
Section 7 demonstrates the qnality of the solutions and Zaremba eigenfnnctions prodnced by 
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the proposed numerical approach, and Section 8, hnally, presents a few concluding remarks. 


2 Preliminaries 

We consider interior and exterior boundary value problems of the form given in equation (1.1) 
for u G (with a Sommerfeld radiation condition in case of exterior problems), where 

n C denotes either a bounded, open and simply-connected domain with a smooth bound¬ 
ary (which we will generically call an “interior” domain) or the complement of the closure 
of such a domain (an “exterior” domain), where the Dirichlet and Neumann boundary por¬ 
tions r D and Tn (h = rDUrAr) are disjoint relatively-open subsets of T of positive measure 
relative to F. Here, calling the ball of radius R centered at the origin, we have denoted 
by 

= {u:ue H\nr] Bn) for all R > O} (2.1) 

the local Sobolev space of order one; of course for bounded sets H. The 

Dirichlet and Neumann data / and g in (1.1), in turn, are elements of certain Sobolev spaces 
(cf. Remark 2.1 item ii) which we dehne in what follows. To do this we follow [23] and we 
hrst dehne, for a given relatively open subset S' C T, the space 

H^/\S) = : supp uCS,ue 

The spaces associated with the Dirichlet and Neumann data / and g are then dehned by 

H^/\S) = : u e 

and, using the prime notation H' to denote the dual space of a given Hilbert space H, 

respectively. 

Remark 2.1. Throughout this paper ... 

(i) ... the term “smooth” is equated to “infinitely differentiable ” and, as indicated above, 
it is assumed that the boundary of the domain D is smooth. 

(ii) .. .it is assumed that the functions f and g in equation (1.1) are smooth. In partic¬ 
ular, it follows that f G and g G — and, thus, that existence and 

uniqueness of the problem (1.1) hold [23]. 

The boundary T can be expressed in the form 

Qn+Qd 

r= U r, ( 2 . 2 ) 

9=1 

where Qd and Qn denote the numbers of smooth connected Dirichlet and Neumann boundary 
portions, and where for 1 < g < Qd (resp. for Qd + 1 < O' < Qd+Qn)i Tg denotes a Dirichlet 
(resp. Neumann) portion of the boundary curve F (see e.g. Figure 1). Clearly, letting 

Jd = {1, • • •, Qd} and .In = {Qd + 1, • • •, Qd + Qn} 
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we have that 


Td = U Tq and = |J Fq 

are the subsets of F upon which Dirichlet and Neumann boundary conditions are enforced, 
respectively. Note that consecutive values of the index q do not necessarily correspond to 
consecutive boundary segments. Throughout this paper it is assumed, however, that no 
Dirichlet-Dirichlet or Neumann-Neumann junctions occur, and, thus, that every endpoint of 
a segment Fq coincides with a Dirichlet-Neumann junction. Clearly this is not a restriction: 
consecutive Dirichlet (resp. Neumann) segments can be combined to produce a partition 
which verihes the assumptions above. 

Remark 2.2. In case D is an exterior domain, problem (1.1) admits unique solutions in 
7F[q^(D). On the other hand, if Q is an interior domain, this problem is not well posed 
for a discrete set of real values kj, j = 1,... ,oo —the squares of which are the Zaremba 
eigenvalues, that is to say, the eigenvalues of the Laplace operator under the correspond¬ 
ing homogeneous mixed Dirichlet-Neumann (Zaremba) boundary conditions (see [23, Th. 

lioj, (Ij). 


3 Boundary integral eqnation formulation 


In what follows we seek solutions of problem (1.1) on the basis of the single-layer potential 
representation 


u = 


Gk{x,y)'ip{y)dsy, 


(3.1) 


where Gk{x,y) = ^H({k\x — y\) is the Helmholtz Green function in two-dimensional space. 
Taking into account well known expressions [9, p. 40] for the jump of the single layer potential 
and its normal derivative across F, the boundary conditions for the exterior (resp. interior) 
boundary value problem ( 1 . 1 ) give rise to the integral equations 


Gk{x,y)ii{y)dsy = f{x) x^Td and 




dGk{x,y) 

dUr 


i){y)dsy = g{x) a; e Fat 


(3.2) 


with 7 = —1 (resp. 7 = 1 ). 

Important properties of both the interior and exterior integral equation problems (3.2) 
relate to existence of eigenvalues of certain interior problems for the Laplace operator under 
either Dirichlet or Zaremba boundary conditions. As shown in what follows, for example. 


1. In case D is an exterior domain the integral equation system (3.2) admits unique 
solutions if and only if k‘^ is not a Dirichlet eigenvalue of the Laplace operator in 

2. For such an exterior domain D the PDF problem (1.1) admits unique solutions for 
any real value of k'^ in spite of the lack of uniqueness implied in point 1 for certain 
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wavenumbers k. A procedure is presented in Appendix A which extends applicability 
of the proposed integral formulation to such values of k. 

3. In case hi is an interior domain, in turn, the integral equation system is uniquely 
solvable provided k"^ is not a Zaremba eigenvalue of the Laplace operator in hi. 

4. The PDE problem (1.1) in such an interior domain hi does not admit unique solutions, 
of course, for values of k for which k'^ is a Zaremba eigenvalue in hi. In this case 
the eigenfunctions of the Zaremba Laplace operator can be expressed in terms of the 
representation formula (3.1) for a certain density 'ijj which satishes (3.2) with / = 0 
and (7 = 0 . 

A detailed treatment concerning points 1, 3 and 4 above is presented in the remainder of 
this section (Theorems 3.3, 3.4 and Dehnition 3.1). A corresponding discussion concerning 
point 2 , in turn, is put forth in Appendix A. 

Definition 3.1. Given an interior (resp. exterior) domain kl and a solution u of (1.1) in 
Hloc{k}), a function w G \ hi) is said to be a solution “conjugate” to u if it satisfies 

Aw + k‘^w = 0 a; e \ , 

I ^ I ^ = r 

w{x) = u[x) a; G 1, 

as well as, in case \ hi is an exterior domain, Sommerfeld’s condition of radiation at 
infinity. Throughout this paper the conjugate solution w in case \ hi is an exterior (resp. 
interior) domain will be denoted by w = Ue (resp. w = Ui). 


Lemma 3.2. The conjugate solutions mentioned in Definition 3.1 exist and are uniquely 
determined in each one of the following two cases: 

1. ^ \ kl is an exterior domain; and 


2. is an interior domain and is not a Dirichlet eigenvalue of the Laplace operator 

in m2 \ n. 


Proof. For both point 1 and point 2 we rely on the fact that the solution u of the problem (1.1) 
is in (see Remark 2.2), and, therefore, by the trace theorem (e.g. [23, Th 3.37]), its 

boundary values he in For point 1 we then invoke [23, Th 9.11] to conclude that a 

uniquely determined conjugate solution w G \ hi) exists, as needed. Point 2 follows 

similarly using [23, Th 4.10] under the assumption that k‘^ is not a Dirichlet eigenvalue in 
the interior domain M^ \ D. □ 


Theorem 3.3. Let D be an exterior domain, and let k eM. be such that kf is not a Dirichlet 
eigenvalue of the Laplace operator in the interior domain M^ \ D. Then the exterior integral 
equation system (3.2) (7 = —1, see also Remark 2.1 item ii) admits a unique solution given 
by 


dui du 
dn dn 


(3.4) 


where u is the solution of the exterior mixed problem (1.1), and where Ui is the (uniquely 
determined) corresponding conjugate solution (Definition 3.1 and Lemma 3.2). 
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Proof. To obtain (3.4) we first consider the Green representation formula for the functions 
Ui and u, 


Ui(x) = I {Gk{x,y)^dsy, a; e \ G, 


u[x) = J \ u- 


ir \ duy 

dGk{x,y) 


- Gk{x,y)^ ] dsy, xeQ, 


(3.5) 


which, in view of the jump relations for the single and double layer potential operators, in 
the limit a; —)■ T leads to the relations 


Ui[X) 


IT \ duy 

u{x) f [ dGk{x,y) 


Gk{x,y)^ ) dsy, a; e r 


dn. 


y 
du 


u- 


dn„ 


(3.6) 


-Gk{x,y )-— dsy, a; e r. 


Since for a; G T^) we have u{x) = Ui{x) = /(a;), the sum of the two equations in (3.6) yields 

= xer„, (3.7) 

and, thus, in view of (3.4), 


f{x)= / Gk{x,y)'4){y)dsy, xeTc,. 


(3.8) 


Similarly, in the limit a; —)■ T the normal derivatives of the integrals in (3.5) give rise to the 
relations 


1 dui 


d 


2 drix drtx Jr 
1 du d f 


^k{x,y)^ ) dsy, a; e r 


driy driy 

JPiPlEl.GAx v)—\ds x€T. 

2 drix dux Jr V dn„ ’ dn„ ' 


i/y !_/ I l/y 

The sum of the equations in (3.9) results in the identity 


1 dui du 


dGk{x,y) ( dui du 


2 drix 2 dux 

or, equivalently, 

du If dui du 


drix 


duy duy ' 


dsy a; G T, 


drix 


2 V dux dria 


du 


+ /f I,er. 

It uTI't 


driy driy j 


(3.9) 


(3.10) 


(3.11) 


But for a; G Tat we have —— = (/(a;), and, thus, equation (3.11) can be made to read 


^(a;) = + f ^^k{x,y) 

2 Jr drix 


f){y)dsy, a;GrAr. 


(3.12) 


7 




Equations (3.8) and (3.12) tell us that the density is a solution of the exterior integral 
equation system (3.2), as claimed. 

In order to establish the solution uniqueness let be a solution of equation (3.2) with 
/ = 0 and g = 0. Since as mentioned above the exterior mixed problem is uniquely solvable, 
the corresponding single layer potential 


V = J Gk{x,y)^{y)dsy (3.13) 

is equal to zero everywhere in hi. It then follows from the continuity of the single layer 
potential that v satishes the Dirichlet problem in the interior domain \ with zero 
boundary values. Since by assumption is not a Dirichlet eigenvalue of the Laplacian in 
it follows that n = 0 in that region as well. Thus, both the interior and exterior normal 
derivatives vanish, and therefore so does their difference The proof is now complete. □ 

Theorem 3.4. Let D be an interior domain. Then we have: 


1 . 


If is not a Zaremba eigenvalue (see Remark 2.2), then the interior integral equation 
system (3.2) ('y=l, see also Remark 2.1 item ii) admits a unique solution, which is 
given by 


^ du dUe 
dn dn 

Here u is the solution of the interior mixed problem (1.1), 
jugate to u (Definition 3.1 and Lemma 3.2). 


(3.14) 

and Up is the solution con- 


2. If k^ is a Zaremba eigenvalue, in turn, any eigenfunction u satisfying (1.1) with / = 0 
and g = 0 can be expressed by means of a single-layer representation 

u{x) = j Gk{x,y) dsy, xeDUT. (3.15) 

where Ug denotes the conjugate solution corresponding to the eigenfunction u (Defini¬ 
tion 3.1 and Lemma 3.2). 


Proof. We hrst consider properties that are common to Zaremba solutions and eigenfunc¬ 
tions, and which therefore relate to both points 1 and 2 in the statement of the theorem. 
For any given solution u of the interior mixed problem (1.1) {u can be either the unique 
solution of the interior mixed problem in the case k/^ is not an eigenvalue, or any eigenfunc¬ 
tion satisfying (1.1) with / = 0 and = 0) the conjugate solution Up is uniquely defined 
(Lemma 3.2). Letting 

w(x) = (3.16) 

using the Green representation formula for u and Up, 


u{x) 

Ue{x) ■■ 


Gk{x,y) 


du dGk{x,y) 


dn,. 


u- 


dn,. 


dGk{x,y) due 

Up -- Gk{x,y) 


dn,. 


dn, 


dSy, 

dSy, 


x G D, 


X eR‘^\Q, 


(3.17) 



and taking into account the jump relations for the double layer potential as well as the fact 
that Ue{x) = u{x) for a; G r we obtain 


u{x) = j Gk{x,y) dsy = w{x), a; G T. (3.18) 

Similarly, taking normal derivatives of both sides of each equation in (3.17) at a point a; G F 
we obtain the equations 


1 du _ d f [^ ( .du dGkix,y)\ 

2 dn^ drix Jr\ ^ ^ duy ^ duy ) 

Idue _ d f f dGk{x,y) .due\ 

2dn^ drix Jr V duy ^ duy) 


dSy^ 


d s 

UjOy, 


a; G r, 
a; G r, 


whose sum yields 


du 

drix 


1 f du 

2 V(9n^ 


<9Me A f dGk{x,y) /^ ^ 

dn^J Jr dn^ \dny duy) ^ dn^ 


(3.19) 


(3.20) 


We now conclude the proof by applying these concepts to points 1 and 2 in the statement 
of the theorem. 


1. In case is not an eigenvalue for the Laplace-Zaremba problem (1.1), equations (3.18) 
and (3.20) evaluated for a; G F^ and a; G Fat, respectively, show that the density tjj 
given by (3.14) satishes the integral equation system (3.2) with 7 = 1, as desired. 

2. In case is an eigenvalue for the Laplace-Zaremba problem (1.1), in turn, let u 
denote a corresponding eigenfunction. Equations (3.18) and (3.20) along with the 
Green representation formula show that u{x) = w{x) for all x E Q. In other words, 
equation (3.15) is satished and the proof follows in this case as well. 

The proof is now complete. □ 


4 Singularities of the solutions of equations ( 1 . 1 ) and ( 3 . 2 ) 

With reference to equation (2.2), let y^ = (|/i,|/ 2 ) G F be a point which separates Dirichlet 
and Neumann regions F^j and Vq^ (gi G Jd and q 2 G Jn) within F. In order to express the 
singular character around y^ of both the solution u{y) of problem (1.1) {y = (|/i,|/ 2 ) G G) 
and the corresponding integral equation density 'iphj) in equation (3.2) {y = (|/i,|/ 2 ) G F) we 
consider the neighborhoods 

G° = Gnfi(|/0,r), F;^ =F,,n 5 ( 1 / 0 , r) and Fj^ = F,, n 5(|/0, r) (4.1) 

of the singular point y^ relative to G, F^j and Tg^, respectively. Here for a set H C A 
denotes the topological closure of H in B{y^, r) denotes the circle centered at y^ of radius 
r, and r > 0 is sufficiently small that B{y^,r) only has nonempty intersections with F^ for 
the Dirichlet index q = qi and the Neumann index q = q2- Additionally, we use certain 
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functions UyO = Uyo[z)^ 'i/'yO = ^lo{d) and -0^0 = ^lo{d) where the Dirichlet (resp. Neumann) 

function (resp ipyo) is the density as a function of the distance d to the point in 

(resp. r°J, and where z = {yi — y^) + i{y 2 — y^) is a complex variable (see Figure 2). The 

functions UyO, %b\ and are given by 
^ y y 

Uyo{z) = u{y), 

^{y) = ^lo{d{y)) 2/erJ^, (4.2) 

V’(2/) =^yo(d(|/)) i/eFj^, 

where, as mentioned above 

Z = {yi - yl) + i{y 2 - yl) and d{y) = - y\Y + (?/2 - ylf. (4.3) 

It is known [29, 30] that, under our assumption that the curve F is globally smooth, for 
any given integer M the function UyO (z) can be expressed in the form 

UyO (z) = \og{z)PyS^ + \og(z)Py^ + Pyo^ + o{z^) (4.4) 

for all z in a neighborhood of the origin, where Pyl^^ Py^ and P^^ are A/'-dependent 
polynomial functions of z, z, z^P, z^P, z\og{z)^ z\og{z). 

Remark 4.1. In the asymptotic expansion (4.4), and, indeed, in all similar asymptotic 
expansions in this paper, it is assumed that none of the right hand side polynomials contain 
terms that, multiplied by the relevant factors, could be included in the error term. 

Under our standing assumption of smoothness of the domain boundary the following two 
theorems provide 1) Finer details on the asymptotics (4.4) as well as 2) A corresponding 
asymptotic expression around y^ for the solutions (d) and (d) of the integral-equation 
system (3.2). 
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Theorem 4.2. Lety^ he a Dirichlet-Neumann point as described above in this section. Then, 
given an arbitrary integer Af, the function UyO (z) can be expressed in the form 

Uyo{z) = + o{z^) (4.5) 

around y^, where is an M-dependent polynomial function of its arguments; see Re¬ 
mark 4-T 

Theorem 4.3. Let y^ be a Dirichlet-Neumann point. Then given an arbitrary integer M the 
functions (d) and {d) can be expressed in the forms 

^lo{d) = + o{d-^~^) and 

around d = 0, where Qy'^ and Qy’o^ are M -dependent polynomial functions of their argu¬ 
ments; see Remark f.l. 

Note in particular that Theorem 4.2 shows that, under our assumptions of boundary 
smoothness, all logarithmic terms in equation (4.4) actually drop out. The proofs of these 
theorems (which are given in Sections 4.2 and 4.3, respectively) utilize a certain conformal 
map introduced in Section 4.1 that transforms 12° into a semicircular region. 


4.1 Conformal mapping 

Following [29], in order to establish Theorem 4.2 we identify with the complex plane 
C via the aforementioned relationship z = {yi — yf) + i{y 2 — y^) ^ {.Hi — 2 / 15 I /2 — 2 / 2)1 
and we utilize a conformal map z = w{f) which maps the semi-circular region = 
G C : 1^1 < 4 and Im(^) < 0} ) in the complex Cplane (Figure 3) onto the domain 12° 
(equation (4.1)) in the complex plane. We assume, as we may, that w maps the origin to 
itself and that the intervals {Im(^) = 0, 0 < Re(^) < A} and {Im(^) = 0, —A < Re(^) < 0} 
are mapped onto the boundary segments F°^ and F^^, respectively. 



Figure 3: Semi-circular and semi-annular Green-identity regions. 


Letting 


R(0 = Uyo{w{f)) 


(4.7) 
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we note that, in view of the relation A^U{^) = AzUyo{w{^)) ■ |w'(0)P satisfied by a complex 
analytic function w (see [15, eq. 5.4.17]), U satishes the second order elliptic equation 


AU + K{^)U = 0 

for 

f e int(T>A), 

(4.8) 

U{0 = F{0 

for 

lm(0 = 0,Re(O > 0 , 

(4,9) 


for 

Im(,^) = 0,Re(.^) < 0, and 

(4.10) 

U(0 = M(0 

for 

lei = A 

(4.11) 


Here F{^) = f{w{^)), G{^) = g{w{^)) and M{^) = u{w{^)). (The function M is thus 
obtained from the restriction of the solution u to the set (911° \ (r°^ blTg^); see equation (4.1) 
and Figure 2). 


4.2 Proof of theorem 4.2 


The proof of this theorem, which, under the present scope of smooth-domain problems 
establishes a result stronger than [29, Th. 3.2], does incorporate some of the lines of the 
proof provided in that reference. In what follows we use the Laplace-Zaremba Green function 




log|t-^| +log|t-^| -2 log 




2 log 




(4.12) 


for the lower half plane with homogeneous Dirichlet (resp. Neumann) boundary conditions 
on the positive (resp.negative) real axis in terms of the complex variables t = ti + it 2 = 
(GjG) and ^ = ('Ci)'C 2 )- The branches of the square roots in (4.12) are given by 

Vi = Vi = V~fV^ ^ where (ptVt) and denote 

polar coordinates in the complex t- and .^-plane respectively (—tt < 9t,9^ < n). Note that, 
with these conventions the domain Da in the t variables is given by pt < 4 and —tt < 9t < 0. 
The following Lemma establishes certain important properties of the aforementioned Green 
function. 


Lemma 4.4. The function H = (equation (4.12 )) is indeed a Laplace-Zaremba Green 

function for the lower half plane with a Dirichlet-Neumann junction at the origin, that is we 
have AtH(t,f) = —S(t — f) and H satisfies 

jATT 

Hit,f) = i} for = 0 and - —(t,.^)=0 for 9t = —7i. (4-13) 

OUt 

In addition, for a certain constant G we have 



Vwt,o) 


dt<G 


(4.14) 


for all f ^ C. 


Proof. The function H{t, f) can be re-expressed in the form 


H{t,0 


1 , \Vi-Vi\\Vi+Vi\ 

— log-^. 

\Vi + ViWVi - ^/i\ 


(4.15) 
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The first statement in equation (4.13) follows from the relations \y/i — a/^| = \ Vi — a/^| and 
Iv^ + v^l = I V^+ V^l, which hold for 6'^ = 0 (t > 0) since, in view of our selection of branch 
cuts we have 

^ = (4.16) 

In order to establish the second statement in (4.13) and equation (4.14) we consider the 
relations 


d 

— loglVf- (2:1 + *2:2)1 = 
^loglv/f - ( 2:1 + *^ 2)1 


_+_ 

2v^—ti {^zl + (v^—+ ^ 2 )^) 


_ ^2 _ 

2 v^ ((^1 — + ^2) 


for 9t = —TT, and 

for 9t = 0, 


(4.17) 


which are valid for all complex numbers z = Zi -\- iz 2 , and we note that, on the axis t 2 = 0 
we have ^ Thus, the second statement in (4.13) follows from application of the first 

equation in (4.17) to each one of the four logarithmic terms that result from expansion of 
equation (4.15). In order to establish a bound of the form (4.14), hnally, we use the second 
equation in (4.17) to obtain the expression 




IMV?)! 


+ 


V(v^ +Re(y^))2 + Im(y^)2 (^_ Re(y^))2 + im(v/^)^ 

(4.18) 

for the absolute value of the derivative of (4.15) with respect to t 2 - It is then easy to check 
that 



7rlm(y^) 


arctan 


- Re( Ve) 
Im(y^) 


+ arctan 


y/A + Re(v^) A 
Im(y^) j 


(4.19) 


for all G C \ [0; A]. Since the right-hand side of equation (4.19) is uniformly bounded for 
G C \ [0; A] and since for ^ G [0; A] the integrand in (4.14) vanishes (in view of the second 
expression in (4.17)), we see that there exists a constant C such that equation (4.14) holds 
for all G C, as needed, and the proof is thus complete. □ 

The proof of the theorem 4.2 is based on a bootstrapping argument which is initiated by 
the simple but suboptimal asymptotic relation put forth in the following lemma 


Lemma 4.5. The solution U of the problem (4.8)-(4.11) satisfies the asymptotic relation 

U{0 = o{e) (4.20) 

for all —4 < /i < 0. 

Proof. To establish this relation we consider the Green formula 
7/(0 = j H{t,OAU{t)dxtdyt + !^U{t)-^H{t,0-H{t,0-^U{t)^dst. (4.21) 
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Since H satisfies (4.13) as it befits a Green function for (4.8)-(4.11), denoting by the 
radius-yf part of ODa it follows that 


1^(01 < 



H{t,OK{t)U{t)dxtdyt 


Da 


+ 


F{x)-^H{t,C)dt 




+ 


Ta 




(4.22) 


For the integral over the outer arc F^ in (4.22) we have 



orit 




<C, 


for 1^1 < A/2, 


(4.23) 


as it can be checked easily in view of the boundedness of the integrands for ^ near the origin. 
Taking into account that u G (see Remark 2.2), on the other hand, it easily follows 

that U G H^{Da), and thus, bounding the absolute value of the hrst integral in (4.22) by 
means of the Cauchy-Schwarz inequality, for ^ near 0 we obtain the uniform estimate 



H{t,OK{t)U{t)dxtdyt 


Da 


<\\H\\l 2 \\U\\l 2 max(iF). 


(4.24) 


In order to estimate the second and third integrals in (4.22) we note that the integrals of 
the absolute values of H and dH/drit are uniformly bounded, as can be checked easily for 
the former, and as it is established in Lemma 4.4 of the latter. The boundedness of F and G 
(which are smooth functions in view of Remark 2.1) thus implies the uniform boundedness 
of the function U near the origin. The relation (4.20) therefore follows for all — | < /i < 0, 
and the proof is complete. □ 


Corollary 4.6. The derivatives of the solution U of the problem (4.8)-(4.11) satisfy the 
asymptotie relation 

D^U = o {C-^) (4.25) 

for all —4 < /i < 0. 

Proof. See [29, Section 4]. □ 

A key element in the bootstrapping algorithm mentioned at the beginning of this section 
is a representation formula for the function U that is presented in the following Lemma. 


Lemma 4.7. The solution U of equation (4.8)-(4.11) admits the representation 


U{0 = -^{Ai {-K{t)U {t) , e, 1) + Ai {-K{t)U {t) , e, 1) - 2 Ai {-K{t)U {t) , 2) 

- 2Ai {-K{t)U (t) ,l2)- A3(F(t), e, 1) - As {F{t),l l) + 2 A 3 {F{t),f, 2) 

+ 2 A 3 {F{t),l 2) - A 2 (G i-t) , -e, 1) - A 2 {G i-t) , -e, 1) + 2 A 2 (G i-t) , -e, 2) 

-I- 2A2 (G (—t), —f, 2)} -I- Pi + P2 ’ 

(4.26) 
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where 


Ai(g(t),^,77) := 


a A 

' 

rA 


g(t)log ptdptdOt, 


Mdit)A,v) ■= / (l{t)\og tv dt, 


(4.27) 


Mq{t)A,v)-=J q{t)--^^^og tv -dt, 

and where pi and p 2 denote power series with positive radii of convergence. 
Proof. Applying the Green formula on the set 

-Da ,(5 = {'C ^ C : (5 < |.^| < A and Im(,^) < 0} 

(right portion of Figure 3) we obtain the expression 


(4.28) 


D(0 = 


H(t,f)AU(t)dxtdyt + 


Uit)—H{t,0-Hit,0^Uit)U.s,. 


Further, for hxed 7^ 0 we have 


H{t, ^) = O (yi^ as f —)■ 0, 


and therefore, in view of Lemma 4.5 and (4.25), 

u{t)H{t,o = o{tn, 


U{t)-yH{t,i) - H{t,f)-yU{t) = o (D-i) 


(4.29) 


(4.30) 


(4.31) 


Letting F5 denote the radius-(5 arc within the boundary ODa^s of Da,s (Figure 3), and noting 
that for t G F^ we have A = A, in view of (4.31) we obtain 


^ iij. V vv v/i vvv- V/uctiix 

U{t) ^H{t, f) - H{t, I dst ^ 0 as <5 ^ 0. 

uTlt C/Tlt J 


(4.32) 


Further, exploiting the fact that the Green function (4.12) is a jointly analytic function of 
\/f and \/f for t G Fa and f around .^ = 0, we obtain 


L ~ =p. (A) +P2 (v^). 

where pi and p 2 denote power series with positive radii of convergence. 
Letting 5 —)■ 0 in (4.29) and using (4.32) and (4.33) we hnally obtain 


(4.33) 


U{0 = 


H{t,Oi-Km{t))ptdptdet- J F{x)-—H{t,0dt 


G{-t)H{-t,f)dt + pi +P2 (^\f^ ■ 


(4.34) 


Using the dehnitions (4.27) for the functions Ai, A2 and A3, equation (4.34) is equivalent to 
equation (4.26) and the proof is complete. □ 
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In order to determine the singular character of U{^) around = 0 (and therefore that 
of UyO (z) around = 0 ) we study the corresponding asymptotic behavior of each one of the 
A-terms in equation (4.26). An important part of this discussion is the following Lemma, 
which presents certain regularity properties of the operators Ai, A 2 and A 3 . 

Lemma 4.8. Let a > 0, (3 > —1 and 7 > —1, and let r] = 1 or rj = 2. Then 

Ai [tn\ e, r,) = + Pi (e") + P 2 (e ) , (4.35) 


A 2 (t^ e, h) = and (4.36) 

A3(f“,e,h) = C^ 6 r + C^7r+P5(e'/’')+P6(f^’'). (4.37) 

For general functions g(t) G C^{Da) and h(t) G C^((0,A]) satisfying g(t) = o{tF) and 
hit) = o(f") as f —)■ 0 , further, we have 

Ai {g{t), e, V) = qi + q 2 (e) + o(e+2), (4.38) 

A 2 (p(t),^, 7 ) = g 3 +g 4 +o(r^^) and (4.39) 

As {h{t),f, ri) = gs j + o{C), (4.40) 


(see Remark 4.1). Here pi (resp. qi), i = 1,... , 6 , are power series with positive radii of 
convergence (resp. polynomials), Ci, i = 1.. .7 are complex constants, and the expansions are 
i-times differentiable as f ^ 0— in the sense of Wigley: the derivatives of the left hand sides 
in (4.38) through (4.40) are equal to the corresponding derivatives of the first two term of 
the right hand sides, with error terms given by the ‘formal” derivatives of the corresponding 
error terms — e.g. d/df 

Proof of Lemma 4-8. The proof follows by specializing the proofs of Lemmas 7.1 — 7.2 and 
8.1- 8.6 in [29]. □ 

We are now ready to provide the main proof of this section. 

Proof of Theorem 4-2. Since the solutions UyO and U are related by equation (4.7), using the 
classical result [27, Th. IV] (which establishes, in particular, that the conformal mapping w is 
smooth up to and including the boundary for any smooth segment of the domain boundary; 
see also [26]) and expanding w{f) in Taylor series around .^ = 0, we see that it suffices to 
prove that for an arbitrary integer Ai we have the asymptotic relation 

(7K) = S"K'yA5+oK") (4.41) 

for the solution Lf{f) of the problem (4.8)-(4.11), where = Q''^(r,s) is a polynomial 
function of the independent variables r and s. 

The proof now proceeds inductively. The induction basis is given by the asymptotics (4.20) 
of the function U. To complete the proof we thus need to establish the following inductive 
step: provided that for some integer C the function U can be expressed in the form 

U{0 = t'^) + o(e^-") as e ^ 0 , (4.42) 
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for some 0 < A < 1, where V‘^ = V'^{r,s) is a polynomial function of the independent 
variables r and s, then a similar relation holds for U with an error of order and 

for a certain polynomial C ): 

U{i) = (4.43) 

To do this we apply Lemma (4.8) to each term on the right hand side of equation (4.26). 
For the terms including the operator Ai, for example, such an asymptotic representation 
with error of the order o(.^'^’''^“^) can be obtained by using the assumption (4.42) and the 
£-th order Taylor expansion of the smooth function K{t) around t = 0, and by applying 
equations (4.35) with = 0,1/2,1,..., £, 7 = 0,1/2,1,..., £ and equation (4.38) with 7 = 
£ — A to the resulting polynomial and error terms for the product K{t)U{t). The terms that 
contain the operators A2 and A3 can be treated similarly on the basis of Taylor expansions 
of the functions F(t) and G(t) around t = 0 and application of equations (4.36), (4.37) with 
/3 = 1,..., £, q; = 1,..., £ , and equations (4.39), (4.40) with 7 = £ — A and a = £ — A + 1. 
The inductive step and therefore the proof of Theorem 4.2 are thus complete. □ 

4.3 Proof of theorem 4.3 

The relationship between the density ip and the PDF solution u is given by equation (3.4) if 
is an exterior domain and by equation (3.14) if is an interior domain. Throughout this 
section we assume that is an interior domain, and, thus, that tp is given by equation (3.14); 
the proof for exterior domains is analogous. 

In order to establish the singular character of the density tp we hrst seek an asymptotic 
expression for the conjugate solution Ue near z = 0 (see Definition 3.1). Using a conformal 
mapping approach for Ue similar to the one used in Section 4.1 for the solution u of the 
problem (1.1), in this case we employ a conformal map = n(.^) which maps the semi-circular 
region Da depicted in Figure 3 in the complex ,^-plane onto the domain B{y^,r) \ in the 
complex z—plane (see Figure 2). We assume, as we may, that v maps the origin to itself 
and that the intervals {Im(,^) = 0,0 < Re(.^) < A} and {Im(,^) = 0, —A < Re(.^) < 0} are 
mapped onto the boundary segments F°^ and F^^, respectively (see equation (4.1)). Following 
Section 4.1 in this case we introduce the function U(.^) = Ue{v{p,)) and we note that V satishes 
the second order elliptic problem (cf. [15, eq. 5.4.17]) 


AV + Ki{^)V = 0 

for G int(iA2i), 

(4.44) 

U(0 = uMO) 

for Im(,^) = 0, and 

(4.45) 

V{0 = Mi{0 

for 1^1 = A, 

(4.46) 


where Mi is given by Mi{p,) = Ue{v{(,)). 

The following Lemma parallels Lemma 4.5. 

Lemma 4.9. The solution V of the problem (4.44)-(4.46) satisfies the asymptotic relation 

U(0=o(n (4.47) 

for all —4 < p < 0. 
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Proof. Employing the Laplace Green function 

Giii, 1^ - ^1 - log |t - ^ I} (4.48) 

for the Dirichlet problem (4.44)-(4.46) and applying the Green formula to the functions V 
and Gi on the domain we obtain 




Gi{t,f)AV{t)dxtdyt + 




(4.49) 

Since Gi(t,f) = 0 for Im(,^) = 0 and since ODa = [—A, A] UL^i, the triangle inequality yields 


11^(01 < 


Gi{t,f)Ki{t)V{t)dxtdyt + ^j V{t)—Gi{t,f)dt 

V(t)AG,(t,()-Gi(t,()Av(t)\ds, . 


(4.50) 


For the integral over the outer arc in (4.50) we have 


I^V{t)Aa,{t,()-G,{t,()Av(t)'^ds, <Ci tor |«|<^/2, (4.51) 

where Gi is a nonnegative constant, as it can be checked easily in view of the boundedness 
of the integrands for ^ near the origin. From Lemma 3.2, further, it easily follows that 
V G H^{Da), and thus, bounding the absolute value of the first integral in equation (4.50) 
by means of the Gauchy-Schwarz inequality we obtain the bound 

[ [ Gi{t,f)Ki{t)V{t)dxtdyt < ||Gi||l 2 (^^)||F||l 2 (o^) max(iLi(t)) (4.52) 

for all ^ G Da- As is well known, hnally, double layer potentials for bounded densities are 
uniformly bounded in all of space (see e.g. [16, Lemma 3.20]). It follows that the second 
integral in equation (4.50) is uniformly bounded for G since, in view of (4.45), Defini¬ 
tion 3.1 and Theorem 4.2, F is a bounded function for t 2 = Im(t) = 0. The relation (4.47) 
thus follows for all — 4 < /i < 0 and the proof is complete. □ 

Corollary 4.10. The derivatives of the solution V of the problem (4.44)-(4.46) satisfy the 
asymptotie relation 

D’^V = o (4.53) 

for all —4 < /i < 0. 

Proof. See [29, Section 4] □ 

We now proceed with the main proof of this section, which is based on an inductive 
argument similar to the one used in the proof of Theorem 4.2. 
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Proof of Theorem f.3. Applying the Green formula on the set Da ,5 (equation (4.28)) and 
letting (5 —)■ 0 we obtain 

^(0= [ [ i-K,it)Vit))G,it,Odt- r Vit)^G,it,0dt + Pii0+P2i0: (4.54) 

Kf Kf D Kf — A. i 

where pi and p 2 denote power series with positive radii of convergence. In view of equa¬ 
tion (4.45), Dehnition 3.1 and Theorem 4.2, on the other hand, we see that, for any given 
integer £, the boundary values oiV at .^2 = 0 satisfy 

V{iiA) = V{0=Ue{0=Vi^{e'\t") + o{e) for e 2 = Im(0 = 0 (4.55) 

(see (4.5)). Relying on equations (4.54) and (4.55) as well as Lemma 4.8, an inductive 
argument similar to the one used in the proof of Theorem 4.2 shows that for any integer M 
the function V satishes an asymptotic relation of the form 

4"(0 = + o(e^) as e ^ 0 , (4.56) 

where V'^ is an A/”—dependent polynomial. In view of corollaries 4.6 and 4.10, substitution 
of the normal derivatives of equations (4.56) and (4.5) for Im(,^) = 0 into equation (3.14) 
yields 

for e 2 = Im(O = 0, (4.57) 

where are A/”—dependent polynomials. The desired asymptotic relations (4.6) now follow 
by re-expressing (4.57) in terms of the distance function d, and the proof is thus complete. □ 

5 Singularity resolution via Fourier Continuation 

Theorem 4.2 tells us that the solutions of Zaremba problem (1.1) possess a very specihc sin¬ 
gularity structure near the Dirichlet-Neumann junctions—which, as shown in Theorem 4.3, 
are inherited by the solutions of the corresponding integral equation system (3.2). In partic¬ 
ular, equation (4.6) shows that the integral equation solutions can be expressed as a product 
of the function l/\/~d and a smooth function of v^, where d denotes the distance to the 
Dirichlet-Neumann junction. 

The question thus arises as to how to incorporate the singular characteristics of the 
integral equation solutions in order to design a numerical integration method of high order 
of accuracy for the numerical discretization of the integral equation system (3.2). A relevant 
reference in these regards is provided by the contribution [5] (see also [31]), which provides 
a high-order solver for the problem of scattering by open arcs. As is known, the open- 
arc integral equation solutions possess singularities around the end-points: they can be 
expressed as a product of the function 1 / y/d and a smooth function of d —or, in other words, 
the asymptotics of the integral solutions are functions which only eontain powers of y/d with 
exponents equal to {2n — 1) for n > 0. In particular [5] a change of variables of the form 
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t = coss, 0 < s < TT in parameter space completely regularizes the problem, and it thus 
gives rise to spectrally accurate numerical approximations of the form 

n 

tjj ~ Cj cos(js) 0 < S < TT (5.1) 

i=o 

for the integral-equation solutions 

As shown in Theorem (4.3), on the other hand, the asymptotic expansions of the integral- 
equation solutions Ip considered in this paper contain all integer powers of '/d, and therefore, 
as established in [ 1 ], a cosine change of variables such as the one considered above leads to 
a full Fourier series—containing all 27r-periodic cosines and sines, 

n 

tp ~ Cj cos{js) + Dj sin(js) 0 < s < tt, (5.2) 

j=0 

even though values for tp can only be determined for 0 < s < tt. The key element that allows 
such expansions in the extended interval [0, 27r] is the Fourier Continuation (FC) method 
introduced in [3,6] and hrst suitably generalized to the present context in [1]. This leads 
to a Fourier series that converges with high-order accuracy to the integral density tp in the 
interval [ 0 , 7 r]. 

Exploiting such rapidly convergent Fourier expansions, a numerical method for Zaremba 
boundary value problems is presented in the following section. 

6 Numerical algorithms 

Closed form expressions are presented in [ 2 ] for the integrals of products trigonometric func¬ 
tions and logarithms that appear in equation (3.2) upon substitution of the expansion (5.2); 
as shown in [ 1 ], use of such expressions leads to highly accurate approximations of the integral 
operators in equation (3.2). In particular, these algorithms rely on Nystrom discretization 
of the solution tp] in view of the structure of the integrand, the discretization points used are 
given as the image of a uniform mesh in s variable under the change of variables t = cos(s). 
The Fourier series (5.2) is obtained via application of the Fourier Continuation method in 
the s variable, and a discrete version of the integral equation system is thus obtained on the 
basis of a resulting matrix A. 

These discrete operators were used in reference [1] as important building blocks of an 
efficient algorithm for evaluation of eigenvalues and eigenfunctions of the Laplace operator 
under Zaremba boundary conditions. Figure 4 presents high-frequency Zaremba eigenfunc¬ 
tions obtained by means of that solver. In what follows we use these discrete operators as 
well as the vector f of values of given functions / and g at the aforementioned discretization 
points to produce solutions of the Zaremba problem (1.1). For simplicity we rely on an 
Lt/-decomposition applied to the linear system Ac = f to obtain a discrete approximation 
c of the solution tp. Note, however, that as mentioned in Theorem 3.3, the integral equation 
system (3.2) is not invertible for a certain discrete set of values of k (spurious resonances) 
that correspond to Dirichlet eigenvalues on the complementary set \ hi. A numerical 
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methodology described in Appendix A extends applicability of the proposed solvers to all 
frequencies, including spurious resonances. A variety of numerical results obtained by means 
of the aforementioned Zaremba boundary-value solvers are presented in the following section. 



Figure 4: High frequency Zaremba eigenfunctions. 


7 Numerical results 

This section presents results of applications of the new solvers to problems of scattering by 
smooth obstacles under boundary conditions of Zaremba type. This entails solution of the 
problem (1.1) for exterior domains H and for which the right hand sides / and g are given 
by 


p _ ^ika-x _ ^2/!:(cos(a)3^i+sin(a)3^2)) 

g = n^- Ve , 

where a is the angle of incidence. 

In our hrst experiment we apply the solver to a kite-shaped domain whose boundary is 
given by the parametrization 

Xi = cos(t)-|-0.65 cos(2t) — 0.65 and 0:2 = 1.5sin(t), (7.2) 

assuming the Neumann and Dirichlet boundary portions Tjv and Tjj corresponding to t G 
[7r/2;37r/2] and its complement, respectively. Figure 5 depicts the scattered and total helds 
that result as a wave with wavenumber k = 40 and incidence angle a = n/8 impinges on 
the object. Figure 6 demonstrates the high-order convergence results for the value of the 
scattered held u{xo) at the point xq = (1, 2 ) in the exterior of the domain. 
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Figure 5: Scattering from a kite-shaped domain under Zaremba boundary conditions. Left: 
Scattered field. Right; Total held. 



Figure 6: Convergence of the value M(a;o) (xq = (1, 2)) for a kite shaped domain with k = 10. 

A similar example concerns scattering by the unit disc under Dirichlet and Neumann 
boundary conditions prescribed on the left and right halves of the disc boundary, respectively. 
Figure 7 displays the scattered and total helds that result as a wave with wavenumber /c = 50 
and incidence angle a = 7i/8 impinges on the disc. 
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Figure 7: Scattering from a disc under Zaremba boundary conditions. Left: Scattered field. 
Right: Total field. 

To conclude this section we present a brief comparison of the proposed solvers with one 
of the most efficient Zaremba solvers previously available [18,19]. This method, which is 
based on iterative inverse preconditioning and which applies to a variety of singular prob¬ 
lems, has been implemented in a numerical MATLAB package which is freely available [20]. 
Unfortunately, Zaremba boundary conditions are not implemented in that package. At any 
rate, numerical experiments suggest that the execution times required by the algorithm [20] 
are not as favorable as those required by the solvers proposed in this paper. Indeed, the 
solver [20] (which is applicable to domains with corners) required a computing time of 0.46 
seconds to solve the Dirichlet problem for Helmholtz equation with wavenumber fc = 2 on 
the unit disc (whose boundary was partitioned into two Dirichlet arcs) by mean of GMRES 
iterations with a GMRES residual of 10“^^, while a GMRES-based implementation of the 
FG-based solver presented in this paper runs in 0.06 seconds for the significantly more chal¬ 
lenging Zaremba problem for the Helmholtz equation on the same domain, with the same 
incident wave frequency and with the same GMRES residual. (These and all other numerical 
results presented in this section were obtained on a single core of a 2.4 GHz Intel E5-2665 
processor.) Such time differences, a factor of eight in this case, can be very signihcant in 
practice, in contexts where thousands or even tens of thousands of solutions are necessary, as 
is the case in inverse problems as well as in our own solution [1] of high-frequency eigenvalue 
problems, etc. The main reason for the difference in execution times is that even though the 
iterative solver [20] requires a limited number of iterations, certain iteration-dependent ma¬ 
trix entries occur in this solver (in view of corresponding iteration dependent discretization 
points it uses), which require large number of evaluations of expensive Hankel functions at 
each iteration, and, thus, a signihcant computing cost per iteration. 

8 Conclusions 

This paper presented, for the hrst time, detailed asymptotic expansions near Dirichlet- 
Neumann junctions for solutions of Zaremba problems on smooth two dimensional domains. 
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By precisely accounting for the singularities of the boundary densities and kernels on the ba¬ 
sis of Fourier-Continuation expansions, fnrther, discretizations of high-order accuracy were 
obtained for relevant boundary integral operators. The resulting integral-eqnation solver 
allows for accurate and efficient approximation of the highly-singular Zaremba solutions at 
both the low- and high-frequency regimes. 
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A Appendix: Exterior solution at interior resonances 

This appendix describes an algorithm for evaluation of the solution of the problem (1.1) 
for an exterior domain hi, and for a valne of k'^ that either eqnals or is close to an interior 
Dirichlet eigenvalue of the Laplace operator in the bonnded set \ hi. As mentioned in 
Section 3 in this case the system of integral equations (3.2) does not a have a nnique solntion. 
However, the solution of the PDF is nniqnely solvable for any valne of k. 

The non-invertibility of the aforementioned continnous systems of integral eqnations at a 
wavennmber k = k* manifests itself at the discrete level in non-invertibility or ill-conditioning 
of the system matrix A := A{k) for valnes of k close to k*. Therefore, for k near k* the 
nnmerical solntion of the Zaremba problems nnder consideration (which in what follows will 
be denoted by u := Uk{x) to make explicit the solntion dependence on the parameter k) 
cannot be obtained via direct solntion the linear system A{k) ■ c = f . As is known, however, 
the solntions u = Uk oi the continnous boundary valne problem are analytic fnnctions of 
k for all real valnes of k —inclnding, in particnlar, for k eqnal to any one of the spnrious 
resonances mentioned above—and therefore, the approximate valnes Uk{x) for k snfficiently 
far from k* can be used, via analytic continnation, to obtain corresponding approximations 
aronnd k = k* and even at a spnrious resonance k = k*. 

In order to implement this strategy for a given valne of k = ko it is necessary for onr 
algorithm to possess a capability to perform two steps: 

1. To determine whether ko is “snfficiently far” from any one of the spnrious resonances k*. 

2a. If ko is “snfficiently far” then simply invert the linear system by means of either an LU 
decomposition or the usnally already available Singnlar Valne Decomposition (which 
is nsed to determine the “distance from resonance”). 

2b. If ko is not “snfficiently far” from one of the spnrions resonances k*, then obtain 
the PDE solntion at ko by analytic continnation from solntions for valnes of /c in a 
neighborhood of ko which are “snfficiently far” from k*. 

Here k is said to be “snfficiently far” from the set of spnrions resonances provided the 
corresponding linear system can be inverted withont signihcant error amplihcations. It has 
been noticed in practice [24] that the regions within which inversion is not possible are 
very small indeed, in snch a way that analytic continnation from “snfficiently far” can be 
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performed to the singular or nearly singular frequency with any desired accuracy. For full 
details in these regards see [24], 

Numerical results conhrming highly accurate evaluation of the PDF solution even for 
resonant frequencies are presented in Figure 8 for the case of the FC-based solver applied 
to the Zaremba boundary value problem on the unit disc. The solution errors are displayed 
for two frequencies: k = 11 (where the solutions are obtained using an LU decomposition) 
and the resonant frequency k = 11.791534439014281 (with solutions obtained by means of 
analytic continuation). Clearly, the proposed approach can tackle the spurious-resonance 
problem without difficulty. 



n, number of boundary di.scrctization points 
Figure 8: Convergence comparison at a regular and a resonant frequency. 
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